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A spatial stochastic model is developed which describes the 3D 
nanomorphology of composite materials, being blends of two different 
(organic and inorganic) solid phases. Such materials are used, for ex- 
ample, in photoactive layers of hybrid polymer zinc oxide solar cells. 
The model is based on ideas from stochastic geometry and spatial 
statistics. Its parameters are fitted to image data gained by elec- 
tron tomography (ET), where adaptive thresholding and stochastic 
segmentation have been used to represent morphological features of 
the considered ET data by unions of overlapping spheres. Their mid- 
points are modeled by a stack of 2D point processes with a suitably 
chosen correlation structure, whereas a moving-average procedure is 
used to add the radii of spheres. The model is validated by comparing 
physically relevant characteristics of real and simulated data, like the 
efficiency of exciton quenching, which is important for the generation 
of charges and their transport toward the electrodes. 

1. Introduction. Using methods from stochastic geometry and spatial 
statistics, a stochastic model is developed which describes the 3D nanomor- 
phology of composite materials, being blends of two different (organic and 
inorganic) solid phases. Such materials are used, for example, in photoac- 
tive layers of hybrid polymer zinc oxide (ZnO) solar cells where the two 
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solid phases play the role of a polymeric electron donor, consisting of, for 
example, poly(3-hexylthiophene), and an inorganic ZnO-electron acceptor, 
respectively. There is a great advantage of polymer solar cells due to their 
potentially low production costs, in comparison with classical silicon solar 
cells. However, the efficiency of polymer solar cells critically depends on the 
intimacy of mixing of the donor and acceptor semiconductors used in these 
devices to create charges as well as on the presence of unhindered percola- 
tion pathways in the individual solid phases of the composite material to 
transport positive and negative charges toward electrodes; see, for example, 
Yang and Loos (2007). It is therefore very important to have tools at one's 
disposal which are suitable to analyze and model the 3D morphology of these 
materials quantitatively. So far, no such tools are available in literature due 
to the fact that imaging of the 3D morphology in high resolution is a difficult 
task. The first 3D images of photoactive layers in polymer solar cells, gained 
by means of electron tomography (ET), have been published only recently; 
see van Bavel et al. (2009) and Oosterhout et al. (2009). 

In the present paper, such 3D images are used to fit our model to real data. 
The model then helps to get a better insight into the impact of morphology 
on the performance of polymer solar cells and, simultaneously, it can be used 
for virtual scenario analyses, where model-based morphologies of solar cells 
are simulated to identify polymer solar cells with improved nanostructures. 

The model developed in this paper is based on methods from stochas- 
tic geometry and spatial statistics; see Kendall and Molchanov (2010) and 
Gelfand et al. (2010) for comprehensive surveys on recent results in these 
fields. In particular, stationary marked point processes are considered as 
models for complex point patterns extracted from ET images, where the 
points are associated with additional information, so-called "marks." 

Note that point processes in 3D have been used for many years to ana- 
lyze geometrically complex point patterns; see, for example, Baddeley et al. 
(1987). More recently, further case studies in 3D point process modeling 
have been performed, for example, in Ballani, Daley and Stoyan (2005), Beil 
et al. (2005) and Stoica, Gregori and Mateu (2005); see also Baddeley et al. 
(2006). Besides, there are many monographs dealing with point processes 
in multidimensional spaces and their statistical inference and simulation. 
We refer, for instance, to Daley and Vere- Jones (2008), Diggle (2003), Illian 
et al. (2008), M0ller and Waagepetersen (2004), as well as Stoyan, Kendall 
and Mecke (1995). 

The paper is organized as follows. Section 2 briefly describes the consid- 
ered solar cells and the corresponding image data on which the model is 
based. In particular, in Section 2.3, the main ideas of a multi-scale approach 
to the segmentation of 3D images are summarized, which has been devel- 
oped recently in Thiedmann et al. (2011). The crucial step of this approach 
is to find an efficient representation of the binarized and morphologically 
smoothed images by unions of overlapping spheres. 
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Then, in Section 3, the spatial stochastic model for the ZnO phase is 
introduced, separately for morphologically smoothed ZnO domains (macro- 
scale) and for those parts representing the difference between the smoothed 
and nonsmoothed binary images (micro-scale). Based on unions of overlap- 
ping spheres representing the ZnO domains, that is, marked point patterns 
extracted from ET images, a stochastic model is built for the smoothed 3D 
morphology (macro-scale) of the photoactive layers considered in this paper. 
Since a strong correlation of midpoints of spheres in z-direction is observed, 
we propose a multi-layer approach considering sequences of correlated 2D 
point processes to model the 3D point patterns of midpoints. The members 
of these sequences belong to a suitably chosen class of planar Poisson clus- 
ter processes, being parallel to the x-y-plane. In particular, a generalized 
version of Matern cluster processes is considered, where the cluster points 
are scattered in uniformly oriented ellipses around their cluster centers (Sec- 
tion 3.1.1). To model the 3D point patterns of midpoints, a Markov chain 
with stationary initial distribution is constructed, which consists of highly 
correlated Matern cluster processes (Section 3.1.2). It can be seen as a sta- 
tionary point process in 3D, where the radii of spheres are considered as 
marks of this point process (Section 3.1.3). Subsequently, a spatial stochas- 
tic model for the micro-scale part of the morphological structure is developed 
(Section 3.2). It is used to invert morphological smoothing and completes 
our model for the 3D morphology of hybrid polymer-ZnO solar cells. Fur- 
thermore, a method to fit model parameters to real image data is proposed. 

Section 4 deals with model validation. To evaluate the goodness of fit, 
we compare model characteristics which have been computed from real and 
simulated data, respectively, like the volume fractions of voxels contribut- 
ing to monotonous percolation pathways through the photoactive layer, the 
distribution of spherical contact distances, and the probabilities of exciton 
quenching. These characteristics have already been considered in Oosterhout 
et al. (2009), since they are strongly related with the performance of solar 
cells. Finally, Section 5 concludes and provides a short outlook regarding 
possible future research. 

2. Polymer solar cells. In this section some necessary background in- 
formation regarding the functionality of polymer solar cells is provided, to- 
gether with corresponding image data on which the model is based. 

2.1. Photoactive layers. We consider photoactive layers of hybrid poly- 
mer zinc oxide (ZnO) solar cells where the two solid phases play the role 
of a polymeric electron donor, consisting of, for example, poly(3-hexylthio- 
phene), and an inorganic ZnO-electron acceptor, respectively. Upon expo- 
sure to light, photons are absorbed in the polymer phase and so-called "ex- 
citons," that is, photoexcited electron-hole pairs, evolve. Excitons are neu- 
tral quasi-particles which diffuse inside the polymer phase within a limited 
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Fig. 1. Schematic layout of a polymer- ZnO thin film solar cell, showing the percolation 
of photogenerated holes (+) and electrons (—). 

lifetime; see Shaw, Ruseckas and Samuel (2008). If an exciton reaches the in- 
terface to the ZnO phase, it is split up into a free electron (negative charge) 
in the ZnO and a hole (positive charge) in the polymer phase. This process 
is commonly referred to as quenching, because it reduces the intrinsic fluo- 
rescent decay of the exciton in the polymer. Provided that the electrons in 
the ZnO phase and the holes in the polymer phase reach the electrodes at 
the top and bottom of the photoactive layer, respectively, photocurrent is 
generated. A schematic illustration of the morphology of photoactive layers 
in hybrid polymer-ZnO solar cells is shown in Figure 1, where the electrodes 
are supposed to be parallel to the x-y-plane. For further information about 
polymer solar cells and the physical processes therein we refer, for example, 
to Brabec, Scherf and Dyakonov (2008). 

Note that the extent of blending of the two materials has a large impact 
on the efficiency of these solar cells, because not all excitons are quenched 
due to their limited lifetimes. Thus, a morphology as displayed in Figure 1, 
where both materials are mixed intimately, is desirable since the excitons are 
likely to reach the interface and charges can be generated. In other words, for 
a morphology which would be ideal with respect to this aspect of function- 
ality, each location of the polymer phase should have a distance to the ZnO 
phase that is smaller than the diffusion length of excitons. For each location 
within the polymer phase, the fraction of excitons reaching the interface is 
called the quenching probability at this location. The mean of these quench- 
ing probabilities, that is, the quenching probability at a randomly chosen 
location of the polymer phase, is called the quenching efficiency. 

Furthermore, the existence of unhindered percolation pathways within 
both phases, ZnO and polymer, is crucial since the generated charges have 
to be transported to the electrodes throughout the phases. Because of the 
electric field between the electrodes, these pathways should be preferably 
monotonous. Hence, to obtain solar cells with high efficiency, an intimately 
mixed morphology with monotonous percolation pathways for both charge 
carriers is desirable and should be taken into account when producing de- 
vices. The stochastic model developed in the present paper will be used to 
identify morphologies with improved efficiency by generating virtual mor- 
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phologies and investigating the transport processes of electrons and excitons, 
respectively. This will be the subject of a forthcoming paper. 

2.2. Electron tomography images. The image data have been gained by 
electron tomography (ET); see van Bavel et al. (2009) and Oosterhout et al. 
(2009). In particular, we consider images for three devices with different 
thicknesses of the photoactive layers: 57 nm, 100 nm and 167 nm. For each 
of the three thicknesses, the 3D ET images are given as stacks of 2D grayscale 
images (being parallel to the x-y-plane, say), which are numbered according 
to their location in z-direction. The sizes of these images in the x-y-plane 
are 934 x 911 voxels for the 57 nm film, and 942 x 911 voxels for the other 
two thicknesses. Each voxel represents a cube with side length of 0.71 nm. 

Figure 2 shows representative 2D slices for the three film thicknesses, 
where the darker parts of the images represent the ZnO phase due to a higher 
electron density compared to polymer. The images displayed in Figure 2 
indicate clear structural differences for the three films. With increasing layer 
thickness, the separated domains of polymer and ZnO are getting finer. In 
particular, the thinnest film, that is, the photoactive layer with thickness of 
57 nm, features large domains of both polymer and ZnO. The stochastic 3D 
model to be fitted takes these morphological differences into account. More 
precisely, the model type will be the same for all three layer thicknesses. 
Only the values of some model parameters will be different for the varying 
morphologies; see Section 3 below. 



^^^^^^^ 













Fig. 2. 2D images of hybrid polymer- ZnO solar cells; first column; 57 nm, second: 
100 nm, third: 167 nm; first row: original grayscale images, second: binarized images. 
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To develop a stochastic model for the 3D morphology of hybrid polymer- 
ZnO solar cells, the 3D ET grayscale images have to be binarized appro- 
priately. Binarization is necessary since we need to decide which voxels are 
classified as polymer and which as ZnO. An elementary approach to binarize 
grayscale images is to use a global threshold: voxels are set to white (poly- 
mer) if their grayscale value exceeds a certain threshold, and are otherwise 
set to black (ZnO). However, it is difficult to find a single global threshold 
to binarize the ET images because of the irregular brightness of these im- 
ages. Thus, instead of considering global thresholds, a method of adaptive 
thresholding has been used for binarization; see Thiedmann et al. (2011). 
This method is based on techniques of Yanowitz and Bruckstein (1989) and 
Blayvas, Bruckstein and Kimmel (2006), where the main idea is to construct 
a threshold surface which is location-dependent and takes local conditions 
like overexposure or underexposure into account. Examples of binarizing the 
ET images by adaptive thresholding are displayed in Figure 2. 

2.3. Segmentation of binarized images. In this section we briefly sum- 
marize the main ideas of a multi-scale approach to the segmentation of 3D 
images, which has been developed recently in Thiedmann et al. (2011). The 
crucial step is to find an efficient representation of the binarized and mor- 
phologically smoothed images by unions of overlapping spheres. 

Let B denote the ZnO phase of the binarized images. Since the morphol- 
ogy of the set B is rather complex (see Figure 2), it is difficult to describe 
this morphology directly, just by a single stochastic model. We therefore 
developed a multi-scale approach to represent the ZnO phase by differ- 
ent structural components. Each of them will be described separately by 
suitably chosen stochastic models. More precisely, we distinguish between 
a macro-scale component of the binarized ET images, which is obtained by 
morphological smoothing, and several micro-scale components, which consist 
of those voxels that have been misspecified by the morphological smooth- 
ing; see Figure 3. The intention of morphological smoothing is to reduce 
the structural complexity of the binarized ET images, that is, to omit very 
fine structural components such as "thin ZnO branches," that is, thin ZnO 
parts connected to larger ZnO domains, "isolated ZnO particles," that is, 
small ZnO particles in the polymer domains, and "polymeric holes," that is, 
small polymeric particles inside the ZnO domains. The morphological trans- 
formations which we use for smoothing the ZnO phase B of the original 
binarized ET images are twofold: "dilation" and "erosion." The morpholog- 
ically smoothed version of the set B will thus be denoted by B". 

In the next step a stochastic algorithm [see Thiedmann et al. (2011)] is 
used to efficiently represent the set B" by a union of spheres, which is de- 
noted by B'" . This leads to an enormous data reduction. Another advantage 
of this representation of the set B" by unions of spheres is that it allows the 
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Fig. 3. Original image split up into structural components at two different scales. 



interpretation of the morphologically smoothed ZnO phase as a realization 
of a marked point process where the midpoints of the spheres are the points 
and the corresponding radii the marks. 

Recall that in the macro-scale component B" of the binarized ET images 
as well as in its representation B'" by unions of spheres, some structural de- 
tails of the original ZnO phase B, like isolated particles, thin branches and 
polymeric holes, are omitted. Furthermore, the boundaries of ZnO domains 
are morphologically smoothed and slightly enlarged by dilation. Hence, when 
comparing the sets B and B'" , we observe that some voxels are misspeci- 
fied, that is, indicated as ZnO although originally being polymer, and vice 
versa. The set BAB'" = {BUB'") \ {BnB'") of misspecified voxels is subdi- 
vided into several subcomponents, where each of these subcomponents will 
be modeled separately. First, two main types of misspecifications are dis- 
tinguished: outer misspecifications and inner misspecifications; see Figure 4. 
Each ZnO voxel that is not covered by a sphere, that is, belonging to the 




Fig. 4. First: binarized (nonsmoothed) 2D slice of 57 nm file; second: representation by 
unions of spheres; third: outer misspecifications; fourth: inner misspecifications (boundary 
and interior). 
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set B \ B'" and therefore constituted as polymer, is called an outer mis- 
specification. Typically, thin branches and isolated ZnO particles are outer 
misspecifications. On the other hand, each polymer voxel which is covered 
by a sphere, that is, belonging to B^CiB'" and constituted as ZnO, is called 
an inner misspecification. Inner misspecifications are further subdivided into 
boundary misspecifications and interior misspecifications. On the one hand, 
polymer voxels (belonging to B^), located near the boundary dB" of the 
macro-scale component B" and covered by a sphere, are called boundary 
misspecifications. On the other hand, each inner misspecification which is 
not a boundary misspecification is called an interior misspecification. Typi- 
cally, polymeric holes belong to interior misspecifications. 

3. Stochastic modeling. We now present our approach to stochastic mod- 
eling of the 3D nanomorphology of the ZnO phase in photoactive layers with 
three different thicknesses which are given by the binarized ET images de- 
scribed in Section 2.2. Note that the model type is the same for all three 
thicknesses, just the fitted values of some model parameters are different. 
This means, in particular, that our model can be used for computer-based 
scenario analyses with the general objective of developing improved materi- 
als and technologies for polymer solar cells. 

In accordance with the multi-scale representation of the ZnO phase which 
has been described in Section 2.3, we will establish stochastic simulation 
models separately for the morphologically smoothed ZnO domains repre- 
sented by unions of overlapping spheres, that is, the macro-scale represen- 
tation of the ZnO phase, and for the three types of misspecifications, that 
is, the micro-scale components. 

3.1. Point-process model for systems of overlapping spheres. To begin 
with, we develop a point-process model which describes the macro-scale 
component of the ZnO phase represented by unions of spheres. This model 
is constructed in several steps. First we consider 2D point processes for those 
midpoints of spheres which belong to single slices of voxels, being parallel to 
the x-y-plane. Since a strong correlation of midpoints spheres in z-direction 
is observed, we propose a multi-layer approach considering sequences of cor- 
related 2D point processes to model the 3D point patterns of midpoints. 

The members of these sequences belong to a suitably chosen class of planar 
Poisson cluster processes. In particular, elliptical Matern cluster processes 
are considered, where the cluster points are scattered in ellipses of uniformly 
distributed orientation around their cluster centers. To model the 3D point 
patterns of midpoints, a Markov chain with stationary initial distribution is 
constructed, which consists of highly correlated Matern cluster processes. It 
can be seen as a stationary point process in 3D, where the radii of spheres 
are considered as marks of this point process. The mark correlation func- 
tions, which have been computed for the radii of spheres extracted from 
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ET images, show strong positive correlations for small distances between 
midpoints. Hence, for a given configuration of midpoints, the radii associ- 
ated with these midpoints are not modeled just by independent marking, 
but a certain moving-average procedure is proposed. Note that our model 
of a stationary marked point process in 3D describing the macro-scale com- 
ponent of the ZnO phase is not isotropic. This is in accordance with the 
nanomorphology observed in the ET images; see Oosterhout et al. (2009). 

3.1.1. Elliptical Matern cluster processes. To get an idea which class of 
point processes is suitable to model the midpoints belonging to the individual 
slices of voxels, we consider the pair correlation function g : (0, oo) — ?■ (0, oo) 
of stationary and isotropic point processes in M?. Note that g{r) is pro- 
portional to the relative frequency of point pairs with distance r > from 
each other; see, for example, Illian et al. (2008). Then, for each of the three 
photoactive layers with thicknesses of 57 nm, 100 nm, and 167 nm, the val- 
ues g{r) of the pair correlation function have been estimated for distances r 
within some interval (0,rmax); see Figure 5. 

Since g{r) > 1 for small r > 0, these estimates clearly indicate clustering 
of points, which can also directly be seen from the point patterns shown in 
Figure 5. The clusters appearing in these point patterns seem to be located in 



■) b] cj 




T t T 



Fig. 5. Top: point patterns of midpoints in 2D slices (a, = 57 nm, b = 100 nm, 
c = 167 nm); bottom: estimated pair correlation functions (d = 57 nm, e = 100 nm, 
f = 167 nm). 
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relatively small (bounded) areas, which corresponds to the fact that g{r) ^ 1 
for sufficiently large r > 0. The shapes of the clusters are not circular, but 
rather elongated. Hence, we propose to consider Matern cluster processes, 
where the cluster points are scattered in ellipses of uniformly distributed 
orientation around their cluster centers. 

This class of (elliptical) Matern cluster processes in can be described by 
a vector of four parameters: (Ac, A^, a, b), where Ac is the intensity of the sta- 
tionary Poisson point process {T„, n > 1} of cluster centers, a and b with a > 
6 > are the semi-axes of (random) ellipses £a,b{Tn,Cn) C centered at the 
points Tn of the Poisson process of cluster centers and rotated around 
Tn by random angles Cn which are independent and uniformly distributed on 
the interval [0, vr), and A^ is the intensity of the stationary Poisson processes 
{Sni,i > 1} of cluster members which are released by the cluster centers 
Tn within the ellipses £a^biTn,Cn)- The Matern cluster process is then de- 
fined as the random point pattern {Sn} given by {Sn} = \J'!^=i{{Sni,i > 

<^£a,b{Tn,Cn)), where the sequences {Cn},{Tn},{Sii},{S2i}, ■ ■ ■ are as- 
sumed to be independent. 

3.1.2. Markov chain of Matern cluster processes. There are strong simila- 
rities between consecutive 2D slices in terms of a high correlation of midpoint 
locations in z-direction as well as approximately equal numbers of points. 
Figure 6 shows a series of such consecutive 2D slices from the 57 nm data set. 

As a consequence, it is not suitable to model the stacks of these 2D point 
patterns by sequences of independent Matern processes. But the (vertical) 
correlation structure visualized in Figure 6 can be taken into account by 
considering a Markov chain of Matern processes. This allows us to model 
small displacements of clusters when passing from slice to slice. Furthermore, 
"births" and "deaths" of clusters in z-direction can also be modeled in this 
way. In other words, we consider a certain class of spatial birth-and-death 
processes with random displacement of points; see, for example, M0ller and 
Waagepetersen (2004). 




Fig. 6. Point patterns of midpoints for successive 2D-slices from the 57 nm data set 
fa = slice 35, b = slice 36, c = slice 37). 
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For each integer z >1, let {Bn\n > 1} be a stationary Poisson point 

process in M? with intensity such that < < Ac, and let {6n\n > 1} 
be an independent and identically distributed (i.i.d.) sequence of Bernoulli 
random variables, which is independent of {Bn^}, where P{6n^ = 1) =p for 

some p G (0, 1). Note that {Bn^} and {(^i^^} will be used in order to model 
"births" and "deaths" of cluster centers, respectively. 

(z) 

For each integer z > 1, an i.i.d. sequence {Dn } of random displacement 
vectors d[^^ , D^^ , . . . with values in ig 

considered, which is independent 

of {Si^^} and We assume that the random vectors d[^\D2^\ . . . are 

uniformly distributed in the set b{o, r") \ b{o, r'), where r' and r" denote the 
size of minimum and maximum displacement, respectively; < r' < r". 

(z) 

Then, a (stationary) Markov chain {{Sn },z>l}oi Matern processes can 
be constructed as follows. For z = 1, let {Sn^} be an elliptical Matern cluster 
process as introduced in Section 3.1.1, that is, {Sn^} = U^i({'S'm,* > 1} H 
£a,b{Tn,Cn))- Assume that the "birth rate" A^ and the "survival probabili- 
ty" p satisfy XcP + A^ = Ac, where Ac is the intensity of the Poisson process 
{Tn} of cluster centers. For z = 2, the Poisson process {Tn } of cluster cen- 
ters is then given by {Tj,^\n> 1} = [j-.^w^^iT^ + of'^} U {Bi^\n > 1}. 

The Poisson process {Tn^^} is constructed in the same way as {tP}, that 
is, 

M'\n>l}= U {T^''>+Df}U{B(^\n>l}, 

(2) (3) 

and so on; see also Figure 7. The Matern processes {Sn },{Sn }, ■ ■ ■ are 

built similarly to the construction of {T^}, {T^}, For example, {Sn^} 

is given by 



a) b) c) d] 




Fig. 7. (a) initial 2D point pattern; (b) displacement of cluster centers, including "birth 
and "death;" (c) 3D point pattern; (d) union of spheres. 
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oo 

u(j{{si^l^>l}n£a,M'\&), 

n=l 

where the sequences {Cn^^}, {Sj^^}, {Sj^^}, . . . are defined in the same way as 
{Cn}, {Sii}, {S2i}, ■ ■ ■ introduced in Section 3.1.1. 

The Markov chain {{Sn^},z > 1} of Matern processes introduced above 
can be seen as a stationary point process in 3D. It possesses seven (free) pa- 
rameters: Ac, Xd, CL, b describing its initial distribution, and p, r', r" describing 
the transitions from step to step, whereas the "birth intensity" of (new) 
cluster centers is given by A^ = Ac(l — p). It turned out that suitable choices 
for r',r" are the values of r' = ^2/2 and r" = 1.5. This means that the uni- 

(z) (z) 

form distribution of the displacement vectors D\ , D2 ,■■ ■ is implemented 
as (discrete) uniform distribution on the 8-neighborhood in the considered 
slice of voxels. Techniques for fitting the remaining five parameters of the 
Markov chain {{Sn^},z > 1} are discussed in Section 3.1.4. 



3.1.3. Modeling the radii of spheres. To get an idea which class of mark 
distributions is suitable to model the radii of spheres, we computed his- 
tograms of radii which have been extracted from the ET images for each of 
the three photoactive layers with thicknesses of 57 nm, 100 nm and 167 nm. 
Recall that in the sphere-putting algorithm mentioned in Section 2.3 we only 
consider spheres with a minimum radius of \/3 voxel sizes. Hence, instead 
of computing histograms for the original radii, say, ri,r2, . . . , we computed 
histograms for correspondingly reduced radii r[,r2,..., where r'^ = rn — \/3- 
It turns out that for all three film thicknesses, gamma distributions can be 
fitted quite nicely to the histograms of reduced radii; see Figure 8. The pa- 
rameters k and 9 of these gamma distributions T{k,9) have been estimated 
using the method of moments; see Table 1. 




Fig. 8. Histograms of reduced radii and fitted gamma distributions (solid lines); 
a = 57 nm, b = 100 nm, c = 167 nm. 
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Table 1 

Parameters for gamma distributions of radii 



Parameter 


57 nm film 


100 nm film 


167 nm film 


k 


1.51 


1.36 


1.26 


e 


1.73 


0.88 


0.93 



The mark correlation function of stationary marked point processes is 
considered, which describes the spatial correlations of pairs of marks, de- 
pending on the distance vector of the corresponding pairs of points; see, for 
example, Illian et al. (2008). For each representation of the three photoac- 
tive layers by unions of overlapping spheres, the values K(r) of this function 
have been estimated for distance vectors of length r. They show strong pos- 
itive correlations for radii corresponding to pairs of midpoints with small 
distances from each other; see Figure 9. 

Thus, for a given configuration {s^f\n, z > 1} of midpoints, the radii 
{Rn\'n,z > 1} associated with these midpoints are not modeled just by 
independent marking, but the following moving-average procedure is pro- 
posed. For some m > 1, let {Rn\n, z > 1} be an i.i.d. sequence of T{k/m, 9)- 
distributed random variables, and let (zi, rei), . . . , (zmj ^m) for each index 
(z, n) denote the indices of the m nearest neighbors Sn^\ • • • , Sn^'^ of Sn^ [in- 
cluding the point s^n^ itself]. Then, the radius R^^^ = \/3 + ^i^/^ -h- • • + ^&'^ 
is assigned to the midpoint Sn \ The reduced radius Rn^ — \/3 obtained 
in this way is r(A;, 0)-distributed. It turned out that for m = 4, the esti- 
mated mark correlation functions computed from real data (i.e., original 
point pattern and original radii) show a good resemblance to their simu- 
lated counterparts (i.e., original point pattern and simulated radii) for all 
three thicknesses of photoactive layers; see Figure 9. 




Fig. 9. Comparison of estimated (solid) and simulated (dotted) mark correlation func- 
tions for arbitrary distance vectors m 3D of length r; a = 57 nm, b = 100 nm, c = 167 nm. 
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3.1.4. Model fitting for midpoints of spheres. The (overall) intensity of 
midpoints of spheres, that is, the intensity A of the stationary point pro- 
cess {Sn\z,n>l} can be easily estimated by A = #{5^^^ : 5n^^ G 

(z) (z) 

where #{5n -Sn G W} is the total number of midpoints in the sampling 
window W CM.^ and \W\ denotes the volume of W. Using this formula, 
the following values have been obtained for the spheres extracted from the 
binarized ET images: A = 1.83 • 10-3 for the 57 nm film, A = 5.29 • lO'^ for 
the 100 nm film, and A = 5.15 • 10"^ for the 167 nm film. 

Note that A = Xc^d\Sa,b\, where \£a,b\ denotes the area of an ellipse with 
semi-axes a and b. Therefore, in order to determine A^, the estimator A^ = 
X{Xc\£^-j^\)~^ can be used, provided that an estimator Ac for the intensity Ac 
of cluster centers as well as estimators a and b for the semi-axes a and b 
are available. Similarly to the estimation of A^ discussed above, the estima- 
tor Ac = Ac(l — p) for the birth rate A^, can be considered, provided that 
estimators Ac and p for Ac and p are given. 

Finally, we derive a so-called minimum-contrast estimator for the vector of 
the remaining four parameters Ac, a, b and p, where we traverse the param- 
eter space of these parameters. This means that for each vertex of a certain 
lattice of parameter vectors {Xc,a,b,p), the 3D point process {Sn\z,n > 1} 
of midpoints described in Sections 3.1.1 and 3.1.2 is simulated in the sam- 
pling window W and gamma-distributed radii are added according to the 
moving average procedure described in Section 3.1.3. Then, structural char- 
acteristics of simulated unions of spheres are compared with corresponding 
structural characteristics of unions of spheres extracted from binarized ET 
images. 

In particular, consider the empirical distribution function p^^^ \ [0,oo) — ?• 
[0, 1] of spherical contact distances from polymer to ZnO, computed for 
the macro-scale component of binarized ET images, and let : 
[0, oo) — )• [0, 1] denote the empirical chord-length distribution functions of the 
ZnO domains in these images along the x-, y- and z-axis, respectively. More 
information on spherical contact distance distributions and chord length dis- 
tributions can be found in Ohser and Miicklich (2000) and in Stoyan, Kendall 
and Mecke (1995). Consider the volume fraction V of ZnO, computed for 
the macro-scale component of binarized ET images, and let V' denote the 
volume fraction of those ZnO voxels of the macro-scale component contribut- 
ing to percolation pathways (monotonous and nonmonotonous) through the 
photoactive layer. 

Let F^^,a ^,p, F^f,^ ,^^, ^l!?a,6,p' Pxl]a,b,P corresponding distribution 

functions and Vx^,a,6,p, ^x^abp volume fractions, respectively, obtained 
from simulated 3D morphologies in dependency of Ac, a, b and p. Then, each 
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Table 2 

Parameters for 3D point processes of midpoints 



Parameter 


57 nm film 


100 nm film 


167 nm film 


Ac 


9.0E-5 


1.25E-3 


l.OOE-3 


a 


45 


22 


24 


b 


15 


6 


10 


P 


0.987 


0.991 


0.977 


Ad 


9.59E-3 


lO.OE-3 


6.83E-3 


A'. 


1.17E-6 


1.17E-5 


2.34E-5 



solution {Xc,a,b,p) of the minimization problem 

(Ac, a, b,p) - argmm(u; \\t - i^Ac,aApll + ^ \c,a,b,p\\ 

Xc,a,b,p 

is called a minimum-contrast estimator for (Ac, a, b,p), where w^^^ ,w^^\w^y\ 
w^^^ > and w^"^^ ,w^^ ^ > are some weights such that w^^^ + w^^^ + w^y^ + 
+ yj{v') = ^ud ||F - F|| = sup^gK \F{t) - F {t)\ deuotcs the Kol- 
mogorov distance of F and F. 

The minimization problem described above is solved numerically, that is, 
only a relatively coarse lattice of parameter vectors {Xc,a,b,p) can be taken 
into account. Table 2 summarizes the results obtained in this way, where we 

put 'W^C'D ^ ^(V) ^ ^(V) ^ ^{x) ^ ^(y) ^ ^(z) = 1/(3 . 4) = 1/12. 

The results shown in Table 2 nicely reflect the main structural differences 
and similarities of the patterns of sphere midpoints for the 57 nm, 100 nm, 
and 167 nm films: The estimated values for Ac, a and b indicate that the 
57 nm film has fewer, but larger clusters of midpoints than the 100 nm and 
167 nm films, whereas the intensity A^ of cluster members is similar for all 
three films; see also Figure 5. The survival probability p is close to 1 and, 
therefore, the birth rate AJ. is much smaller than the intensity Ac of "old" 
cluster centers; see Figure 6. 

We also remark that the method of statistical model-fitting explained in 
this section leads to a relatively high degree of visual coincidence between 
simulated and real ET images; see Figure 10. A more formal approach to 
model validation will be given later on in Section 4. 

3.2. Stochastic modeling of clusters of misspecified voxels. We now de- 
velop a modeling approach for the micro-scale part of the morphological 
structure. It is used to stochastically invert the morphological smoothing 
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Fig. 10. Left: 2D slice of the morphologically smoothed 57 nm film, right: 2D slice of 
a simulated union of overlapping spheres, drawn from the fitted 3D model. 



and completes our model for the 3D morphology of hybrid polymer-ZnO 
solar cells. The micro-scale morphology is modeled separately for each of 
the three types of misspecifications, that is, the micro-scale components of 
outer, boundary, and interior misspecifications mentioned in Section 2.3. 

3.2.1. Outer misspecifications. Recall that each ZnO voxel that is not 
covered by a sphere, and therefore constituted as polymer, is said to be an 
outer misspecification. They typically form thin branches or small isolated 
ZnO particles. In the present section, a stochastic model is proposed for 
the locations and sizes of clusters of outer misspecifications, that is, the 
connected components of the set B \ B'" introduced in Section 2.3. We first 
consider this kind of misspecification, because it influences the "correction" 
of boundary misspecifications which will be described in Section 3.2.2 below. 

We assume that the centers of gravity of clusters of outer misspecifications 
form a Cox point process, which is also called a doubly stochastic Poisson 
process in literature. The cluster sizes are considered as marks. In partic- 
ular, under the condition that a realization {(s^f^r^f^)} of the (marked) 
point-process model {{Sn ,Rn ),'n,z > 1} introduced in Section 3.1 is given 
which describes the macro-scale component of the ZnO phase represented 
by the union of spheres ^ = |Jn z>i b{^n\rn^), we assume that the centers of 
gravity of clusters of outer misspecifications can be described by an inhomo- 
geneous Poisson process. Its (conditional) intensity A(x) at location x G 
depends on the distance 6^{x) = mf{\x — y\-y £ C} between x and the union 
of spheres ^, where we put X{x) = if 5^{x) = 0. 

The intensity X{x) at locations x G with 6^{x) > can be estimated 
by analyzing the centers of gravity of clusters of outer misspecifications in 
the set B \ B'" extracted from binarized ET images. In particular, for any 
di,du> with di < du, the (average) intensity \di,du) centers of gravity 
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Table 3 

Model parameters for outer mis specifications 







57 nm 


100 nm 


167 nm 


Intensity 


-^[0,2) 


1.78E- 


3 


5.25E-3 


5.70E-3 




A[2,4) 


5.22E- 


4 


4.31E-4 


6.88E-4 




-^[4,6) 


1.81E- 


4 


1.80E-4 


4.48E-4 




-^[6,8) 


1.04E- 


4 


l.lOE-4 


3.28E-4 




-^[8,10) 


6.90E- 


5 


6.06E-5 


2.31E-4 


Slope 


Q 


-0.90 




-0.67 


-1.13 


Axis intercept 




86.45 




33.62 


30.92 


Variance 




l,889.f 




114.6 


63.9 



at locations x G with b^{x) G \d[^du) can be estimated by 

number of centers of gravity with distance to B'" between di and du 
number of voxels with distance to B'" between di and d^ 

Examples of results for \\di4^) computed from ET images are given in Ta- 
ble 3. For all three film thicknesses the estimated intensity y^\di4i) decreases 
with increasing distance to the macro-scale component B'" of the ZnO phase. 

We assume that within shells around the set ^ = Un z>i ^('^"^^ '""^^) with 
distance to ^ in the distance class [d;, d^), the intensity of centers of gravity of 
outer misspecifications is constant and given by \\di,du) - assume that the 
clusters of outer misspecifications are spheres, the radii of which are given 
in the following way. It turned out that not only the intensity of centers of 
gravity, but also the cluster sizes observed in binarized ET images, depend 
on the distances of centers of gravity to the set B'" . In particular, the cluster 
sizes seem to have a tendency to decrease with increasing distance to the 
set B'" \ see Figure 11, where the estimated mean values for the radii of the 
considered distance classes are shown. 



■I ^ e) 




□EttflFKV In nm Ur»l«EH.v Iri nm Dlftuice in nfn 

Fig. 11. Mean cluster sizes of outer misspecifications, depending on their distance from 
the set B'" ; a = 57 nm, b = 100 nm, c — 167 nm. 
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To integrate this dependency into our simulation model, we fit regres- 
sion lines to the point clouds shown in Figure 11, that is, we assume that 
the points in this figure can be seen as realizations of random variables Yi 
satisfying the linear relation = axi + /3 + e^, where Xi = + d^)/2 is 

the midpoint of the iih distance class [d^f^\d^) and ei is a random error 
term. The parameters a and /3 of this regression line are estimated by the 
method of least squares. As can be seen in the plots of Figure 11, a linear 
model shows just the trend but is not a perfectly fitting model for describing 
the cluster sizes in dependence on their distances to the macro-scale compo- 
nent B'" in the ZnO phase. Hence, we additionally consider the residua £i 
in the linear regression model, where we assume that they follow a normal 
distribution with expectation and variance cr^. The estimated values ob- 
tained for slope a, intercept /3 of y-axis, and variance cr^ of the residua are 
given in Table 3. To ensure a positive size of each simulated cluster, we reject 
negative sizes and generate new realizations as long as a positive cluster size 
is sampled. For simulated clusters of outer misspecifications with a greater 
distance from the set ^ than the intercept of the fitted regression line with 
the X-axis, we put the cluster size equal to zero. This is in accordance with 
real data, because in the binarized ET images such clusters of outer misspec- 
ifications do not occur. In the following, by ^' we will denote a realization 
of the model with included outer misspecifications. 

3.2.2. Boundary misspecifications. After adding the outer misspecifica- 
tions to our model as described in the previous section, we now develop 
an algorithm to remove the so-called boundary misspecifications, which pri- 
marily result from the dilation of the ZnO domains; see Section 2.3. Recall 
that we defined the boundary misspecifications as those misspecified voxels 
within some outer shells of the set B'" , that is, the union of spheres repre- 
senting the morphologically smoothed ZnO phase. The first shell is defined 
as the set of voxels of B'" , with a distance to W\B"' smaller or equal than 1. 
In general, the {i + l)th shell, i = 1, 2, 3, . . . , is defined as the set of voxels 
of B'" with a distance ioW\ B'" m{i,i + l\. 

Table 4 displays the percentage of misspecified boundary voxels in the 
different outer shells of B'" . Note that boundary misspecifications only occur 
in the first outer shell (for the 100 nm and 167 nm films) and in the first 
three outer shells (for the 57 nm film), respectively. As a consequence, the 
simulation model should remove about the same percentage of boundary 
voxels in the corresponding outer shells of ^. 

Some parts of the outer shells of B'" belong to thin branches of ZnO, 
therefore not the complete shells of ^ have to be removed. To include such 
thin branches into the model, we combine the model for the outer misspeci- 
fications introduced in Section 3.2.1 with the following algorithm to remove 
the boundary misspecifications. 
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Table 4 

Fractions of boundary misspecifications in consecutive shells (given in %) 





57 


nm 


100 nm 


167 


nm 


ET data 


Simulation 


ET data 


Simulation 


ET data 


Simulation 


1st shell 


87 


72 


60 


67 


60 


64 


2nd shell 


74 


70 














3rd shell 


71 


68 














4th shell 





















For a given realization ^, we iteratively remove those parts of the outer 
shells of ^ which are not connected to the set ^' \ ^ introduced in the previous 
section. In more detail, to correct the first outer shell, we first determine the 
set of all voxels ?/ C ^ belonging to the first outer shell. Subsequently, all vox- 
els of the set ??i C r/ that are not touching the set ^' that is, whose distance 
to is greater than 1, are removed. Now, the first outer shell is corrected. 
Those parts ?72 = \ ??i of the first outer shell that have not been removed 
since they were located near an outer misspecification are — for technical 
reasons — added to the set \ ^ of simulated outer misspecifications. Hence, 
when correcting the second outer shell, the voxels near the set \ ^) U 772 
are not removed. This reproduces the thin branches as observed in the bi- 
narized ET data. To correct the third shell, the same procedure is repeated. 
The result after additionally adding the boundary misspecifications into the 
model is denoted by ^" . 

3.2.3. Interior misspecifications. As mentioned in Section 2.3, the re- 
maining misspecified voxels in B'" are classified as interior misspecifications. 
These interior misspecifications typically form small polymeric holes inside 
the ZnO domains. 

Our modeling approach for the interior misspecifications is based on the 
assumption that the polymeric holes in the ZnO domains possess spherical 
shapes and are not overlapping, that is, we consider them as hard spheres. 

Similarly to the modeling of the outer misspecifications, we assume that 
the centers of gravity of the interior misspecification clusters form a dou- 
bly stochastic point process, where again the cluster sizes are considered as 
marks. We assume the points of this point process to have a certain min- 
imum distance to each other because the interior misspecifications are 
seen as nonoverlapping spheres. In particular, given the realization ^" , that 
is, a realization of the marked point process for the macro-scale component 
of the ZnO phase introduced in Section 3.1 with included outer and bound- 
ary misspecifications as described in Sections 3.2.1 and 3.2.2, the centers of 
gravity of interior misspecification clusters are assumed to form a (condi- 
tional) Matern hard-core process in ^H^". We assume that the marks of this 
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Table 5 

Estimated parameters for interior misspecifications 



Parameter 


57 nm film 


100 nm film 


167 nm film 


Radius r 


2.50 


1.30 


1.07 


Intensity 


1.37E-3 


5.17E-3 


5.12E-3 



point process are spheres with a constant radius r = that is, the interior 
misspecifications are modeled by nonoverlapping spheres with equal radii. 

The Matern hard-core process in with intensity and hard-core ra- 
dius rfi is a thinned homogeneous Poisson point process, where the remaining 
points have a distance of at least r/j to each other. Further details can be 
found, for example, in Illian et al. (2008). Given the set ^^n^", the centers of 
gravity of interior misspecifications are then modeled by those points of the 
Matern hard-core process which belong to ^ H Hence, this model for the 
interior misspecifications can be called a doubly stochastic Matern hard-core 
process. 

To fit this model to real data, we first estimate the mean volume V of 
the clusters for all three film thicknesses and transform this (mean) volume 
into a radius f of a ball with the same volume V. The corresponding radii 
obtained for the three considered data sets can be seen in Table 5. The 
hard-core radius f/j of the Matern hard-core process is then computed as 
f/j = 2f, which ensures that the spheres of radius f centered at the points 
of the doubly stochastic Matern hard-core process are not overlapping. For 
the intensity of the Matern hard-core process we consider the following 
natural estimator, 

number of disjoint clusters of interior misspecifications 
''^ \B"eb{o,r)\ ' 

where \B" Q b{o,r)\ denotes the volume of the set B" Q b{o,r). The val- 
ues of \h obtained for the three data sets are shown in Table 5. From the 
results given in Table 5 it can be seen that, also with respect to interior 
misspecifications, the 57 nm film behaves rather different than the 100 nm, 
and 167 nm films: in the first case, there are fewer, but larger clusters of 
interior misspecifications than in the latter case. In the following, by we 
denote a realization of our simulation model after including all three types 
of misspecifications. 

Figure 12 shows a realization of the model for the micro-structure, where 
the micro-structure model is applied to the real data, more precisely, the 
representation by a union of spheres B'" of the 57 nm film. The two im- 
ages on the left and the right sides of Figure 12 possess a high degree of 
visual resemblance. See also Section 4 for a more formal approach to model 
validation. 
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Fig. 12. Left: 2D slice of 57 nm film, center: corresponding representation by a union of 
spheres, right: simulated correction using the model for the misspecified voxels applied to 
the representation by a union of spheres. 

4. Model validation. To evaluate the goodness of fit, we compare model 
characteristics which have been computed from real and simulated data, 
respectively. On the one hand, we consider structural characteristics of the 
ZnO nanomorphology like the volume fraction of ZnO, the volume fraction of 
ZnO contributing to monotonous percolation pathways, and the distribution 
of spherical contact distances from polymer to ZnO. On the other hand, we 
consider a physical characteristic, the so-called exciton quenching probabil- 
ity. This characteristic describes the probability that a photo-excited particle 
generates charges. These characteristics have also been used in Oosterhout 
et al. (2009) to characterize the morphology of a photoactive layer, since 
they are closely related to the performance of solar cells. To compare the 
values of these characteristics, obtained from simulated and real data, we 
binarize the ET images using two extreme global thresholds as suggested in 
Oosterhout et al. (2009). Recall that the two global thresholds have been 
chosen in such a way that the ZnO phase can be assumed to be a subset of 
the union of foreground voxels (high threshold) and, vice versa, the polymer 
phase is contained in the union of background voxels (low threshold) . 

It turns out that the estimated values obtained for most of the consid- 
ered image characteristics of these two binarizations can be seen as lower 
and upper bounds, respectively, for corresponding values obtained for sim- 
ulated images. In addition to this, we mention that, in accordance with the 
visual resemblance of images obtained from real and simulated data for the 
macro-scale component (see Figure 10) and the micro-scale component (see 
Figure 12) of the ZnO nanomorphology, the optical resemblance between 
binarized ET images obtained by adaptive thresholding and realizations of 
the complete simulation model is also quite well; see Figure 13. 

4.1. Checking morphological characteristics. For a quantitative valida- 
tion of the stochastic simulation model, we first consider structural char- 
acteristics of the ZnO nanomorphology. For this purpose, we generate 100 
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Fig. 13. 3D cutouts (400 x 400 x 37 voxels) of bmanzed ET images obtained by adaptive 
thresholding (1st row) and realizations of the complete model (2nd row); left: 57 nm film, 
center: 100 nm film, right: 167 nm film. 

realizations of our model, estimate the considered characteristics for each of 
these realizations and compute their mean values. In the case of the spherical 
contact distance distribution function, the pointwise means are considered. 

First, the volume fraction of ZnO is considered, which is one of the most 
important characteristics in structural modeling. The results given in Table 6 
show that for all three film thicknesses, the volume fractions of ZnO com- 
puted from simulated data are between the corresponding bounds obtained 
from the globally thresholded ET images. 

In the next step, the connectivity of the ZnO phase is considered. This 
also is an important characteristic, because only if there is a high connectiv- 
ity, that is, if many percolation pathways exist, the produced charges can be 
transported to the electrodes, where current can be gripped. For estimating 
the connected and monotonously connected volume fractions of ZnO we ap- 
plied the same methods as in Oosterhout et al. (2009). However, in general, 
the (conditional) connected and monotonously connected volume fractions 
of the foreground phase in globally thresholded images do not monotonously 
depend on the values of global thresholds. But as shown in Figure 14, this 
is the case for the globally thresholded ET data of photoactive layers of 
polymer-ZnO solar cells. Hence, also with respect to connected volume frac- 
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Table 6 

Volume fractions of ZnO for globally thresholded and simulated images 



Volume fraction 

Volume fraction Volume fraction with monotonous 







of ZnO 


with connection 


connection 


57 nm 


Low threshold 


0.098 


0.934 


0.872 




Simulated data 


0.112 


0.905 


0.864 




High threshold 


0.172 


0.974 


0.947 


100 nm 


Low threshold 


0.133 


0.890 


0.673 




Simulated data 


0.215 


0.971 


0.910 




High threshold 


0.295 


0.991 


0.936 


167 nm 


Low threshold 


0.128 


0.851 


0.630 




Simulated data 


0.210 


0.943 


0.806 




High threshold 


0.293 


0.979 


0.907 



tions, the values for the two extreme thresholds can be seen as upper and 
lower bounds. With the exception of the 57 nm film, the values for the con- 
nected volume fractions computed from simulated data, given in Table 6, are 
nicely between the corresponding values obtained from globally thresholded 
ET images. Also, the relative errors between the values obtained for the 
adaptively thresholded ET images and simulated images are rather small; 
see Table 7. 

Finally, the spherical contact distribution function F^^^ : [0, oo) — )• [0, 1] 
of the ZnO phase is considered, where F^^^{t) can be interpreted as a con- 
ditional probability that the minimum distance from a randomly chosen 
location to the ZnO phase is smaller or equal than t > 0, provided that the 
considered location belongs to the polymer phase. 
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Fig. 14. Volume fractions of connected (solid lines) and monotonously connected (dashed 
lines) foreground phase, in dependence of the global threshold; a — 57 nm, b = 100 nm, 
c — 167 nm. 
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Table 7 

Volume fractions of ZnO for adaptively thresholded and simulated images 



Volume fraction 

Volume fraction Volume fraction with monotonous 









of ZnO 


with connection 


connection 


57 


nm 


Adaptive threshold 


0.133 


0.963 


0.928 






Simulated data 


0.112 


0.905 


0.864 






Relative error 


-0.158 


-0.060 


-0.069 


100 


nm 


Adaptive threshold 


0.211 


0.980 


0.888 






Simulated data 


0.215 


0.971 


0.910 






Relative error 


0.019 


-0.009 


0.025 


167 


nm 


Adaptive threshold 


0.209 


0.970 


0.851 






Simulated data 


0.210 


0.943 


0.806 






Relative error 


0.005 


-0.028 


-0.053 



Similarly to the situation which we observed for the connected volume 
fractions considered above, it turns out that the spherical contact distri- 
bution functions of the foreground phase in globally thresholded ET im- 
ages depend monotonously on the value of the global threshold. Hence, 
the estimated contact distribution functions Ff^^ and F^^^ obtained for 
the two extreme thresholds can be seen as upper and lower bounds, re- 
spectively; see Figure 15. In addition to this, the spherical contact distri- 
bution F^^^ obtained from simulated data is shown in Figure 15, where 
Fi^'^^{t) < F^'^^{t) < F^'^^{t) for ah considered t > and for ah three layer 
thicknesses. 

In summary, we can conclude that our model fits very well to real data 
regarding the considered structural characteristics. As the model has been 
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Fig. 15. Spherical contact distribution functions. The lower and upper bounds Fi 
and -F^*^^ obtained from globally thresholded ET images are plotted as dashed lines, the 
corresponding results from simulated image data as solid lines; a = 57 nm, b = 100 nm, 
c = 167 nm. 
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developed for analyzing the influence of morphology on the performance of 
solar cells, we also consider a physical characteristic for model validation 
which is described in the following section. 

4.2. Checking probabilities of exciton quenching. Quenching efficiency r/g 
is the probability of a random exciton being quenched; see Section 2.1. It 
is an elementary but important physical characteristic for the efficiency of 
solar cells. 

In a hybrid polymer-ZnO solar cell, absorption of light by the poly- 
mer phase does not directly yield free charge carriers. Instead excitons are 
formed. It is only at the interface of the polymer and ZnO phase that free 
charges are generated by quenching (splitting) of excitons. It is therefore of 
the utmost importance that excitons are able to reach this interface. The ex- 
citon diffusion length in conjugated polymers is typically a few nanometers, 
which puts a considerable constraint on the morphology of polymer solar 
cells. In other words, the efficiency of exciton quenching is very sensitive to 
morphology, making it a suitable way to validate our model. 

Suppose that the polymer phase B'^ = W \ B is given in a cubic sampling 
window C M'^. Then, the overall efficiency rjq of exciton quenching can 
be obtained from the field {n{x),x S B^} of local exciton densities in the 
polymer phase. The exciton density field {n(x),x G B'^} can be computed 
by solving the steady-state diffusion equation [see Oosterhout et al. (2009)] 



where D is the diffusion constant, r is the exciton life time, and g is the 
rate of exciton generation. As a boundary condition we require that all 
excitons at the polymer-ZnO interface be quenched, that is, 77,(3;) = for all 
X G dB^ \ dW. For x G dB^ n dW, cyclic boundary conditions are applied 
in all directions. The exciton lifetime and exciton diffusion rate in P3HT 
are taken from the literature: r = 400 ps and D = 1.8 • 10~^ m^ s~^; see 
Shaw, Ruseckas and Samuel (2008). The rate of exciton generation g is just 
a scaling factor, where we use a value of g = 10^^ m~^ s~^ which is typical 
for 1-sun conditions. 

The diffusion equation is solved numerically to a relative error of less 
than 10"^. Figure 16 shows local exciton density fields {n{x),x £ B'^} for 
adaptively thresholded and simulated data, respectively. 

Once {n(x),x £ B^} is known, the quenching efficiency tjq follows from 
rjQ = 1 — n/{Tg), where n is the average exciton density in the polymer 
domain B'^. Figure 17 compares the quenching efficiencies for original and 
simulated data. The quenching efficiency is also monotonously depending 
on the global threshold. Hence, the values of tjq obtained for the two ex- 
treme thresholds can be seen as lower and upper bounds. The values of the 
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Fig. 16. The local exciton density n{x) normalized to rg for adaptively thresholded ET 
(left) and simulated (right) data. The scalebar, specifying the density of excitons, applies 
to both images. 
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Fig. 17. Quenching efficiencies for globally thresholded and simulated images. 

quenching efficiency for tlie simulated data lie well within these lower and 
upper bounds; see Figure 17. The small relative errors displayed in Table 8 
show that our model fits very well to real data, where quenching efficiencies 
of adaptively thresholded ET images are compared with those of simulated 
data. 



Table 8 

Quenching efficiencies for adaptively thresholded and simulated 

images 





57 nm film 


100 nm film 


167 nm film 


Adaptive threshold 


0.418 


0.794 


0.819 


Simulated data 


0.416 


0.805 


0.834 


Relative error 


-0.010 


0.014 


0.018 
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5. Conclusions and outlook. In the present paper we developed a param- 
eterized stochastic simulation model for the nanostructure of photoactive 
layers of hybrid polymer-ZnO solar cells. The model is based on tools from 
stochastic geometry. Additional to the model itself, we developed a method 
to fit its parameters to real 3D ET image data. 

To establish our model, the adaptively thresholded ET images are seg- 
mented using a stochastic algorithm which consists of two main steps. First, 
the images are morphologically smoothed in order to slightly decrease their 
structural complexity. Then, the morphologically smoothed binary images 
are represented by a system of overlapping spheres, which can be interpreted 
as a realization of a 3D marked point process, where the sphere centers are 
the locations of points and the corresponding radii are their marks. For 
the stochastic simulation model, we use a correlated vector of 2D ellipti- 
cal Matern cluster processes, where the points are subsequently marked to 
create a 3D marked point process. To complete the model, that is, to in- 
clude the structural details which were omitted due to the morphological 
smoothing, a stochastic simulation model for this "micro-scale" component 
is developed afterward. 

As our stochastic simulation model is fully parameterized, we also devel- 
oped techniques for the estimation of the model parameters of all model 
components. Thus, we are able to fit the simulation model to the ET image 
data described in Section 2.2. 

Finally, we validated the simulation model by comparing structural and 
physical characteristics computed from simulated image data with the corre- 
sponding characteristics obtained from globally and adaptively thresholded 
ET images, respectively. In particular, the quenching efficiencies computed 
for realizations of the simulation model agree very well with those of the ET 
images. Hence, our model nicely reflects the diffusion of excitons. 

Since we were able to fit our model to 3D ET data, it has already proved 
its capability to represent realistic nanostructures of photoactive layers of 
hybrid polymer-ZnO solar cells. The fact that the model is parameter-based 
enables us to predict morphologies for film thicknesses, which have not (yet) 
been imaged by 3D ET, by interpolating or extrapolating the fitted model 
parameters. Due to a strong correlation between morphology and efficiency 
of polymer solar cells [see Oosterhout et al. (2009)], the developed simulation 
model is of significant importance for further investigations of polymer solar 
cells. In a forthcoming paper we will also investigate the transport processes 
of charges to the electrodes as described in Koster (2010), additional to the 
structural and physical characteristics considered in the present paper. By 
generating virtual morphologies, which are generated as realizations of the 
developed model with different parameter configurations, and investigating 
the transport processes of electrons and excitons therein, the spatial stochas- 
tic model will be used to identify morphologies of improved efficiency with 
respect to the considered physical characteristics. 
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We also remark that the modehng approach developed in the present 
paper can be applied to various other kinds of image data, including geo- 
graphical data considered, for example, in ecology. Then, for instance, the 
Markov chain of Matern cluster processes introduced in Section 3.1.2 may be 
viewed as a sequence of dependent 2D point processes through time, which 
can model, for example, the temporal movements of species in a given region. 
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